function [ResEuler,yGrid] = EulerEqErrorPruning_demo(model,orderApp,setup)
% Some indices
ny      = model.ny;        %Number of output variables
nx      = model.nx;        %Number of state variables
derivs  = model.derivs;    %Deviatives as stored by Levintal

%Unfolding parameters
names = fieldnames(model.params);
for i=1:length(names)
    eval([names{i},'=',num2str(model.params.(names{i}),16),';']);
end

% The raw grid for the states
xGrid_f = setup.xSim_decom.xt_f;
if orderApp >= 2
    xGrid_s = setup.xSim_decom.xt_s;
end
if orderApp >= 3
    xGrid_rd = setup.xSim_decom.xt_rd;
end
if orderApp >= 4
    xGrid_4th = setup.xSim_decom.xt_4th;
end
if orderApp >= 5
    xGrid_5th = setup.xSim_decom.xt_5th;
end

% The grid for numerical itegration
tmp_e        = cell(1,model.ne);
tmp_w        = cell(1,model.ne);
for i=1:model.ne
    tmp_e(i) = mat2cell(setup.GH_e,length(setup.GH_e),1);
    tmp_w(i) = mat2cell(setup.GH_w,length(setup.GH_w),1);
end
eGrid = cartprod(tmp_e);
wGrid = prod(cartprod(tmp_w),2);
%% Building the residuals
f        = zeros(ny+nx,1);
ResEuler = zeros(size(xGrid_f,2),ny+nx);
yGrid    = zeros(size(xGrid_f,2),ny);

% Setting the dimensions for the state components
xf_cu    = [zeros(nx,1);1];
xs_cu    = [zeros(nx,1);1];
xrd_cu   = [zeros(nx,1);1];
x4th_cu  = [zeros(nx,1);1];
x5th_cu  = [zeros(nx,1);1];

xf_cup   = [zeros(nx,1);1];
xs_cup   = [zeros(nx,1);1];
xrd_cup  = [zeros(nx,1);1];
x4th_cup = [zeros(nx,1);1];
x5th_cup = [zeros(nx,1);1];
for k = 1:size(xGrid_f,2)
    xf_cu    = xGrid_f(:,k);
    if orderApp>=2
        xs_cu  = xGrid_s(:,k);
    end
    if orderApp>=3
        xrd_cu  = xGrid_rd(:,k);
    end
    if orderApp>=4
        x4th_cu  = xGrid_4th(:,k);
    end
    if orderApp>=5
        x5th_cu  = xGrid_5th(:,k);
    end
    % Getting the states at time t
    k_cu  = model.h0(1,1) + xf_cu(1,1) + xs_cu(1,1) + xrd_cu(1,1) + x4th_cu(1,1) + x5th_cu(1,1);
    c_ba1 = model.h0(2,1) + xf_cu(2,1) + xs_cu(2,1) + xrd_cu(2,1) + x4th_cu(2,1) + x5th_cu(2,1);
    a_cu  = model.h0(3,1) + xf_cu(3,1) + xs_cu(3,1) + xrd_cu(3,1) + x4th_cu(3,1) + x5th_cu(3,1);
    d_cu  = model.h0(4,1) + xf_cu(4,1) + xs_cu(4,1) + xrd_cu(4,1) + x4th_cu(4,1) + x5th_cu(4,1);

    % Getting the states at time t+1
    xf_cup(1:end-1,1)     = derivs.hx*xf_cu;
    if orderApp>=2
        xf2_cu            = kron(xf_cu,xf_cu);
        xs_cup(1:end-1,1) = derivs.hx*xs_cu + derivs.hxx*xf2_cu/2;
    end
    if orderApp>=3
        xf3_cu             = kron(xf2_cu,xf_cu);
        xf_xs_cu           = kron(xf_cu,xs_cu);
        xrd_cup(1:end-1,1) = derivs.hx*xrd_cu + derivs.hxx*(2*xf_xs_cu)/2 ...
                            + derivs.hxxx*xf3_cu/6;
    end
    if orderApp>=4
        xf4_cu             = kron(xf3_cu,xf_cu);
        xf2_xs_cu          = kron(xf_cu,xf_xs_cu);
        xs2_cu             = kron(xs_cu,xs_cu);
        xf_xrd_cu          = kron(xf_cu,xrd_cu);
        x4th_cup(1:end-1,1)= derivs.hx*x4th_cu + derivs.hxx*(2*xf_xrd_cu + xs2_cu)/2 ...
                            + derivs.hxxx*(3*xf2_xs_cu)/6 ...
                            + derivs.hxxxx*xf4_cu/24;
    end
    if orderApp>=5
        xf5_cu             = kron(xf4_cu,xf_cu);
        xf3_xs_cu          = kron(xf_cu,xf2_xs_cu);
        xf_xs2_cu          = kron(xf_cu,xs2_cu);
        xf2_xrd_cu         = kron(xf_cu,xf_xrd_cu);
        xs_xrd_cu          = kron(xs_cu,xrd_cu);
        xf_x4th_cu         = kron(xf_cu,x4th_cu);
        x5th_cup(1:end-1,1)= derivs.hx*x5th_cu ...
                            + derivs.hxx*(2*xf_x4th_cu+2*xs_xrd_cu)/2 ...
                            + derivs.hxxx*(3*xf2_xrd_cu+3*xf_xs2_cu)/6 ...
                            + derivs.hxxxx*(4*xf3_xs_cu)/24 ...
                            + derivs.hxxxxx*xf5_cu/120;
    end

    % Getting the controls at time t
    if orderApp==1
        yVec_cu = model.g0 + derivs.gx*xf_cu;
    elseif orderApp==2
        yVec_cu = model.g0 + derivs.gx*(xf_cu+xs_cu)+derivs.gxx*(xf2_cu)/2;
    elseif orderApp==3
        yVec_cu = model.g0 + derivs.gx*(xf_cu+xs_cu+xrd_cu)+derivs.gxx*(xf2_cu+2*xf_xs_cu)/2 ...
            + derivs.gxxx*(xf3_cu)/6;
    elseif orderApp==4
        yVec_cu = model.g0 + derivs.gx*(xf_cu+xs_cu+xrd_cu+x4th_cu) ...
            + derivs.gxx*(xf2_cu+2*xf_xs_cu+2*xf_xrd_cu+xs2_cu)/2 ...
            + derivs.gxxx*(xf3_cu+3*xf2_xs_cu)/6 ...
            + derivs.gxxxx*xf4_cu/24;
    elseif orderApp==5
        yVec_cu = model.g0 + derivs.gx*(xf_cu+xs_cu+xrd_cu+x4th_cu+x5th_cu) ...
            + derivs.gxx*(xf2_cu+2*xf_xs_cu+2*xf_xrd_cu+2*xf_x4th_cu+xs2_cu+2*xs_xrd_cu)/2 ...
            + derivs.gxxx*(xf3_cu+3*xf2_xs_cu+3*xf2_xrd_cu+3*xf_xs2_cu)/6 ...
            + derivs.gxxxx*(xf4_cu+4*xf3_xs_cu)/24 ...
            + derivs.gxxxxx*xf5_cu/120;
    end
    yGrid(k,:) = (yVec_cu-model.g0)';
    % The controls
    c_cu  = yVec_cu(1,1);
    iv_cu = yVec_cu(2,1);
    y_cu  = yVec_cu(3,1);
    la_cu = yVec_cu(4,1);
    n_cu  = yVec_cu(5,1);
    rk_cu = yVec_cu(6,1);
    w_cu  = yVec_cu(7,1);
    
    Euler   = zeros(ny+nx,1);
    xf_cup0 = xf_cup;
    for i=1:length(wGrid)
        % Adding shocks to the states
        xf_cup(1:nx,1) = xf_cup0(1:nx,1) + model.eta*eGrid(i,:)';
        
        % Getting controls at time t+1
        if orderApp==1
            yVec_cup  = model.g0 + derivs.gx*xf_cup;
        elseif orderApp==2
            xf2_cup   = kron(xf_cup,xf_cup);
            yVec_cup  = model.g0 + derivs.gx*(xf_cup+xs_cup)+derivs.gxx*(xf2_cup)/2;
        elseif orderApp==3
            xf2_cup   = kron(xf_cup,xf_cup);
            xf3_cup   = kron(xf2_cup,xf_cup);
            xf_xs_cup = kron(xf_cup,xs_cup);
            yVec_cup  = model.g0 + derivs.gx*(xf_cup+xs_cup+xrd_cup)+derivs.gxx*(xf2_cup+2*xf_xs_cup)/2 ...
                + derivs.gxxx*(xf3_cup)/6;
        elseif orderApp==4
            xf2_cup    = kron(xf_cup,xf_cup);
            xf3_cup    = kron(xf2_cup,xf_cup);
            xf_xs_cup  = kron(xf_cup,xs_cup);
            xf4_cup    = kron(xf3_cup,xf_cup);
            xf2_xs_cup = kron(xf_cup,xf_xs_cup);
            xs2_cup    = kron(xs_cup,xs_cup);
            xf_xrd_cup = kron(xf_cup,xrd_cup);
            yVec_cup   = model.g0 + derivs.gx*(xf_cup+xs_cup+xrd_cup+x4th_cup) ...
                + derivs.gxx*(xf2_cup+2*xf_xs_cup+2*xf_xrd_cup+xs2_cup)/2 ...
                + derivs.gxxx*(xf3_cup+3*xf2_xs_cup)/6 ...
                + derivs.gxxxx*xf4_cup/24;
        elseif orderApp==5
            xf2_cup    = kron(xf_cup,xf_cup);
            xf3_cup    = kron(xf2_cup,xf_cup);
            xf_xs_cup  = kron(xf_cup,xs_cup);
            xf4_cup    = kron(xf3_cup,xf_cup);
            xf2_xs_cup = kron(xf_cup,xf_xs_cup);
            xs2_cup    = kron(xs_cup,xs_cup);
            xf_xrd_cup = kron(xf_cup,xrd_cup);
            xf5_cup    = kron(xf4_cup,xf_cup);
            xf3_xs_cup = kron(xf_cup,xf2_xs_cup);
            xf_xs2_cup = kron(xf_cup,xs2_cup);
            xf2_xrd_cup= kron(xf_cup,xf_xrd_cup);
            xs_xrd_cup = kron(xs_cup,xrd_cup);
            xf_x4th_cup= kron(xf_cup,x4th_cup);
            yVec_cup = model.g0 + derivs.gx*(xf_cup+xs_cup+xrd_cup+x4th_cup+x5th_cup) ...
                + derivs.gxx*(xf2_cup+2*xf_xs_cup+2*xf_xrd_cup+2*xf_x4th_cup+xs2_cup+2*xs_xrd_cup)/2 ...
                + derivs.gxxx*(xf3_cup+3*xf2_xs_cup+3*xf2_xrd_cup+3*xf_xs2_cup)/6 ...
                + derivs.gxxxx*(xf4_cup+4*xf3_xs_cup)/24 ...
                + derivs.gxxxxx*xf5_cup/120;
        end
        c_cup  = yVec_cup(1,1);
        iv_cup = yVec_cup(2,1);
        y_cup  = yVec_cup(3,1);
        la_cup = yVec_cup(4,1);
        n_cup  = yVec_cup(5,1);
        rk_cup = yVec_cup(6,1);
        w_cup  = yVec_cup(7,1);
        
        % Getting the states at time t+1
        k_cup  = model.h0(1,1) + xf_cup(1,1) + xs_cup(1,1) + xrd_cup(1,1) + x4th_cup(1,1) + x5th_cup(1,1);
        c_ba1p = model.h0(2,1) + xf_cup(2,1) + xs_cup(2,1) + xrd_cup(2,1) + x4th_cup(2,1) + x5th_cup(2,1);
        a_cup  = model.h0(3,1) + xf_cup(3,1) + xs_cup(3,1) + xrd_cup(3,1) + x4th_cup(3,1) + x5th_cup(3,1);
        d_cup  = model.h0(4,1) + xf_cup(4,1) + xs_cup(4,1) + xrd_cup(4,1) + x4th_cup(4,1) + x5th_cup(4,1);
        
        %% Model Specific
        % START DISPLAYING f
        f(1,1)=exp(-la_cu)*exp(d_cu)*(1/(exp(c_cu) - B*exp(c_ba1))^ETAc - (B*BETTA*exp(d_cup))/(exp(c_cup) - B*exp(c_cu))^ETAc) - 1;
        f(2,1)=1 - (THETA*exp(-la_cu)*exp(-w_cu)*exp(d_cu))/(1 - exp(n_cu))^ETAl;
        f(3,1)=BETTA*exp(-la_cu)*exp(la_cup)*(exp(rk_cup) - DELTA + 1) - 1;
        f(4,1)=exp(rk_cu) + (exp(a_cu)*exp(n_cu)^ALFA*(ALFA - 1))/exp(k_cu)^ALFA;
        f(5,1)=1 - ALFA*exp(-w_cu)*exp(a_cu)*exp(k_cu)^(1 - ALFA)*exp(n_cu)^(ALFA - 1);
        f(6,1)=1 - exp(-y_cu)*(exp(c_cu) + exp(iv_cu));
        f(7,1)=exp(-y_cu)*exp(a_cu)*exp(k_cu)^(1 - ALFA)*exp(n_cu)^ALFA - 1;
        f(8,1)=1 - exp(-iv_cu)*(exp(k_cup) + exp(k_cu)*(DELTA - 1));
        f(9,1)=1 - exp(-c_cu)*exp(c_ba1p);
        f(10,1)=RHOA*log(exp(a_cu)) - log(exp(a_cup));
        f(11,1)=RHOD*log(exp(d_cu)) - log(exp(d_cup));
        % END DISPLAYING f
        Euler(:,1) = Euler(:,1) + wGrid(i,1)*f(:,1);
    end
    
    ResEuler(k,:) = Euler(:,1)';
end

end